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Abstract 

We report first-principles frozen-phonon calculations for the determination 
of the force-free geometry and the dynamical matrix of the five Raman-active 
A\ g modes in YBa2Cu307. To establish the shape of the phonon potentials 
atomic forces are calculated within the LAPW method. Two different schemes 
- the local density approximation (LDA) and a generalized gradient approx- 
imation (GGA) - are employed for the treatment of electronic exchange and 
correlation effects. We find that in the case of LDA the resulting phonon fre- 
quencies show a deviation from experimental values of approximately -10%. 
Invoking GGA the frequency values are significantly improved and also the 
eigenvectors are in very good agreement with experimental findings. 
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I. INTRODUCTION 



The phonons in high temperature superconductors have been intensively studied in the 
last decade, where Raman-scattering is the most often used experimental technique. The 
phonon frequencies of different symmetries are measured for a variety of high-T c (HTC) 
compounds. On the other hand, much less is known about the corresponding eigenvectors. 
Experimentally, information is accessible when certain atoms in the unit cell are replaced 
by chemically equivalent elements. By analyzing the frequency shifts from Raman measure- 
ments performed on isotope-substituted samples the composition of the normal modes can 
be obtained. The experimentally best investigated HTC crystal is YBa2Cu307_ x . Raman 
measurements on isotope-substituted samples of this material have given insight into the 
coupling of atomic vibrations especially for the two K\ g phonons containing Barium and 
copper oscillations.0! 

In the past years, ab-initio calculations of frequencies and eigenvectors of high-symmetry 
phonons in YBa2Cu307 have been published.H0 These works were carried out by the frozen- 
phonon technique within all-electron band-structure calculations utilizing the local density 
approximation (LDA) for the treatment of electron exchange and correlations. The phonon 
potentials were derived from total-energy values for different ionic configurations. The so 
obtained A lff phonon frequencies show a tendency of being slightly too soft, lying approxi- 
mately 10% below their experimental counterparts. 

This deviance can be discussed in several terms. First, the sole employment of total- 
energy values for establishing the dynamical matrix requires self-consistent electronic- 
structure calculations for a huge number of atomic configurations. When only harmonic 
terms are taken into account this number is proportional to N 2 when N atoms are involved 
in the vibrations of the same symmetry class. Thus the number of actually used total energy 
values hardly exceeds the minimum number of configurations. Anharmonic contributions 
to the phonon potentials which might change the diagonal and off-diagonal elements of the 
dynamical matrix could not be accounted for. Second, the strong electron correlations in 
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high-T c materials are known to be underestimated by LDA leading to some shortcomings in 
the description of the electronic density distribution!. This fact might also bring up certain 
inaccuracies in the description of structural and vibronic properties of this material. 

In this paper we present frozen-phonon calculations employing the linearized augmented 
plane-wave (LAPW) method.i@ We evaluate the atomic forces acting on the ions when 
distorted from equilibrium.!!! Using these data instead of total energy values the geometry 
is optimized and the phonon potentials of the five Ax g modes are constructed. We checked 
the necessity of including anharmonic terms in the description of the energy surface and the 
number of configurations required for a stable fit. This procedure was carried out within 
the local-density approximation first. Alterations of the phonon potentials when treating 
electron exchange and correlations via a generalized gradient approximation (GGA)0 rather 
than via LDA are also investigated. 



II. CRYSTAL STRUCTURE AND SYMMETRY 

YBa2Cu30y is an orthorhombic crystal with the space group Pmmm. Yttrium and the 
chain atoms Cu(l) and 0(1) occupy sites with the full symmetry of the cell (mmm). For 
barium, the copper and oxygen atoms in the planes - Cu(2), 0(2), and 0(3) - and the 
bridging oxygen 0(4) the site symmetry is lower (mm), with a vibrational degree of freedom 
along the c-axis. Thus the five A\ g Raman-active modes are coupled c-axis vibrations of 
these atoms. 



III. METHOD 

A. Energy Surface 

Within the frozen-phonon approach potentials underlying lattice vibrations are tradi- 
tionally constructed through polynomial fitting of total-energy values calculated at selected 
frozen-in distortions of the crystal. If N denotes the number of coupled degrees of freedom 
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in a certain phonon symmetry, the minimal number of total-energy calculations necessary 
for the harmonic approximation of the potentials scales with N 2 . This scaling behavior 
in conjunction with the formidable computational effort of self-consistent band-structure 
calculations for polyatomic structures limits the feasibility of this approach when applied to 
YBa 2 Cu 3 07. Thus, in past publications,i@ the number of actually calculated energy values 
hardly exceeded the minimal number necessary for the quadratic energy surface. In this way, 
neither anharmonic effects in the phonon potentials can fully be taken into account nor is it 
possible to perform a reliable check on the stability of the fitting parameters. These limita- 
tions can be overcome by evaluating the forces acting on atoms when they are displaced from 
equilibrium. Atomic-force formalisms!!!-!] have proved to give an accurate description of the 
total energy gradient with respect to ionic positions and thus enlarge the information on the 
shape of the energy surface obtainable from each frozen-phonon calculation. In fact, for the 
harmonic approximation it suffices to consider a number of distortion patterns scaling with 
N rather than N 2 . The additional numerical effort required for the evaluation of atomic 
forces is not significant. 



B. Atomic Forces 

If the electronic density distribution p(r) in a crystal is precisely known, the Hellmann- 
Feynman theoremBEl states that the atomic force F£* (defined as the negative derivative of 
the total energy E to t with respect to the nuclear coordinate R Q ) is exactly described by the 
electrostatic force exerted on the nucleus a by all other charges of the system: 

t _ dEtat _ v (Rq - R/j) r (R Q - r) 3 

F « = " dR7 " ~%^| R a - R/3 | 3 qa Jn p{r) | R Q - r / T {3A) 

with qp denoting the nuclear charge at site Rg. This electrostatic force is called Hellman- 
Feynman force and is denoted F^ F in the following. In practical ab-initio calculations 
charge densities are never exact and so correction forces have to be added to the HF-force 
to obtain an accurate description of the total-energy gradient^. The calculation of these 
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forces requires an analytic treatment of the first-order change of the energy expression when 
the ionic positions R Q are displaced by a small amount 5H a . The energy expressions depend 
on the approximations made when simplifying the many-body crystal Hamiltonian. In the 
case of the all-electron band-structure method LAPW,§ which is based on density-functional 
theory (DFT) and a local exchange-correlation potential, this results in two corrective force- 
terms so that the atomic force can be written as 

F a J = F* F + F c ° re + F 1 * s (3.2) 

Y c ° re is due to the simplifications made for the core elctrons' wave-functions. The core-states 
are assumed to be non-dispersive in k-space and are obtained by solving the Schrodinger 
equation neglecting the non-spherical part of the effective crystal potential V e ff(r). It was 
shown in Ref. || that ~F c ° re can explicitly be calculated by 



//Ce(OW e //(r)d 3 r (3.3) 



* a = " / P. 

Jn 

F^ 55 is related to the fact that only a finite number of basis functions can be taken into ac- 
count when solving the eigenvalue problem for the valence electrons by a variational scheme. 
In our particular case, the wave-functions are linear combinations of linearized augmented 
plane waves (LAPWs). Since these LAPWs depend on the atomic positions the incomplete- 
ness of the set of basis functions (IBS) leads to one part of the correction force F ! a BS . A 
further contribution termed T) a stems from the discontinuity of the second derivative of the 
wave-functions at the atomic sphere boundaries. The total F IBS is given byi 



K BS = -J2 



W; 



'At + V eff - eM + mf + V ef f - eA 



(3.4) 



with 



D« = ^2 w i ( f k*(r)f^(r)| MTQ - ^(r)f^(r)| Jnt ldA a . (3.5) 

i JMT a 1 J 

ijji denotes the wavefunction of a valence-state, w\ the respective occupation number, and 
Ei the corresponding eigenvalue. The integration in the term D Q runs over the surface of 



the atomic sphere around the nucleus a. ipi(r)\MT a is the eigenstate expressed in terms of 
the local basis functions within the atomic sphere. ipi(r)\i nt is the interstitial part of the 
wave-function where the basis set consists of plane-waves. 

IV. CALCULATIONS 

A. Energy Surface and Dynamical Matrix 

The zone-center A\ g Raman-active modes in YB^CusOy are coupled c-axis vibrations 
of Ba, Cu(2), 0(2), 0(3), and 0(4). Using experimentally determined lattice constants^ 
(a = 3.8231 A, b = 3.8864 A, c = 11.6807 A) atomic forces acting on these atoms were 
calculated for 19 different distortion patterns in the vicinity of the force-free geometry. 10 
patterns corresponded to single-atom displacements, in the remaining nine configurations 
two atoms were distorted simultaneously. The maximum amplitude of displacement from 
equilibrium was approximately 0.05 A. This yielded a set of 95 values for establishing a 
least square fit of the energy surface with respect to the five coupled degrees of freedom. To 
optimize this fit anharmonic coefficients supplemented the linear and quadratic terms where 
necessary. On the basis of the 15 independent harmonic force-constants of the resulting 
polynomial function we set up the dynamical matrix for calculating phonon frequencies 
and eigenvectors. Additionally, the stability of the fit parameters was tested by repeating 
the fitting procedure on the sole basis of the 10 single-atom displacement patterns. The 
resulting change in the phonon frequencies was very small (less than ±2 cm -1 ) indicating a 
high reliability of the polynomial coefficients. 

B. Computational Details 

The self-consistent band-structure calculations for each frozen-in distortion were done 
using the full-potential LAPW method as implemented in the WIEN95-code.i Valence 
states and semicore states (Y-4s, Y-4p, Ba-5s, Cu-3p) were treated in the same energy 
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window by using local orbitalsEj as an extension to the LAPW basis set. The expansion of 
the wave functions included 1950 LAPW's (RK max = 7.1). The sphere radii were chosen to 
be 2.5 a.u. for Y and Ba, 1.8 a.u. for Cu(l), 1.9 a.u. for Cu(2) and 1.55 a.u. for 0(1), 0(2), 
0(3) and 0(4). Charge densities and potentials in the atomic spheres were represented by 
spherical harmonics up to L = 6, whereas in the interstitial region these quantities were 
expanded in a Fourier series with 3500 stars of K in the LDA case and 5000 stars of K 
in the GGA case. The radial functions of each LAPW were calculated up to I = 12 and 
the non-spherical potential contribution to the Hamilton matrix (Gaunt coefficients) had 
an upper limit of / = 6. For Brillouin-zone (BZ) integrations a 12 x 12 x 4 k-point mesh 
was used yielding 72 k-points in the irreducible wedge. 

C. Correction Terms 

The unreliability of the pure HF-force for calculating the atomic force within the LAPW 
method has been demonstrated by Yu et alJ for phonons in simple systems as Mo and Si. In 
the case of YBa2Cus07 we can confirm the importance of the corrective force terms for an 
accurate description of the total-energy gradient. In Table | we list the various contributions 
to the atomic forces when the two atoms 0(2) and 0(4) are displaced simultaneously from 
equilibrium. With the exception of 0(2), F IBS and F core are of the same magnitude as F HF 
but of different signs. Due to this cancellation the resulting total force ( |3.2|) is one or even 
two orders of magnitude smaller than its constituents. Thus a very precise calculation is 
required in order to obtain reliable results. 

V. RESULTS 
A. LDA Results 

For the case, where the electronic exchange and correlation effects are treated by LDA, 
the energy surface is well described by harmonic terms plus a cubic term in the 0(4) co- 
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ordinate. The minimum energy positions of this polynomial function are zs a = 0.1812 c, 
z Cu(2) — 0.3530 c, Zo(2) = 0.3789 c, Zom = 0.3783 c, and Zqu\ = 0.1586 c, which 
agree well with the experimental structure parameters (zbu = 0.1843 c, z Cu (2) = 0.3556 c, 
z O{2) = 0.3772 c, Zo(3) = 0.3789 c, and zo(4) = 0.1584 c). The eigenvalues and eigen- 
vectors of the corresponding dynamical matrix are displayed in Table [TJ. The two low-lying 
modes are governed by hybridized Ba-Cu motions. Further we observe an appreciable ad- 
mixture of 0(2), 0(3) and 0(4) oscillations in the 330, 387 and 452 cm -1 modes. The 
importance of the anharmonic term in the 0(4) coordinate is reflected in the fact, that its 
exclusion from the fitting procedure leads to a lowering of the corresponding apex mode (452 
cm" 1 ) by 8 cmT x . The introduction of other diagonal or off-diagonal anharmonic terms had 
the maximum effect of ±1 cmT 1 on the phonon frequencies. 

In the comparison of the calculated phonon-frequencies with their experimental counter- 
parts it can be noticed that the discrepancies lie in the range of —10%. The only exception 
is the out-of-phase motion of 0(2) and 0(3) where the agreement with the experimental 
frequency is excellent. These findings are close to the results of Rodriguez et al.0 based on 
LMTO total-energy calculations but reveal larger differences to LAPW frozen-phonon cal- 
culations by Cohen et al.i particularly for the three modes dominated by the motions of the 
oxygens. They find a frequency of 513 cm" 1 for the apex mode, 12% above our calculated 
value. For the 0(2)-0(3) mode and the 0(2)+0(3) mode their values are 312 cm -1 and 
361 cm -1 , whereas our results are 335 cm -1 and 387 cm -1 , respectively. As little informa- 
tion is given on the computational details in the paper of Cohen et al.i we can not give a 
definit explanation for the discrepancies. We can state, that the influence of anharmonicities 
neither on the diagonal nor on the off-diagonal elements of the dynamical matrix is strong 
enough to account for the different values. Parts of the differences might be connected to the 
fact, that Cohen et. al used experimental equilibrium positions to optimize each coordinate 
separately. It should be noted that in a recent refinement!^ of their work the frequency of 
the apex oxygen mode is in much better agreement with our results. 
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B. GGA Results and Discussion 



It is well known that the correlation effects between the electrons in Y-Ba-Cu-0 com- 
pounds seem to play an important role for the electronic properties of these materials: Mea- 
surements of the electric-field gradient in YBa2Cu307 have revealed that band-structure 
calculations on the basis of LDA show certain inaccuracies in the description of the ma- 
terials' electronic charge distribution.! In addition, such calculations fail completely in re- 
producing the insulating and antiferromagnetic ground state of YBa2Cu3O6-0 This raises 
the question if correlation effects also influence the structural properties of YBa 2 Cu 3 7 and 
thus might account for the theoretical underestimation of the Ax g phonon frequencies. To 
obtain an insight into the shortcomings of LDA concerning vibronic properties, we used a 
generalized gradient approximation proposed by Perdew et al.0 for the treatment of the 
exchange-correlation energy and potential which is now calculated from the local electronic 
density as well as from its gradient. We established the energy surface of the A lg phonons 
following the same procedure as in the LDA case. 

We find that the GGA scheme enhances the diagonal terms in the dynamical matrix by 
approximately 15% (Ba), 10% (Cu), 3% (0(2)), 5% (0(3)), and 4% (0(4)), respectively in 
comparison to the LDA results. On the other hand, the coupling between the copper and 
the three oxygen atoms is significantly decreased. The same applies to the force-constants 
coupling the apex vibration to the motion of the other two oxygens. The GGA approach 
neither changes the anharmonic cubic term in the 0(4) coordinate noticeably, nor does 
it cause further anharmonicities in the energy surface. Resulting phonon frequencies and 
eigenvectors are displayed in Table |ITl[ We observe a significant improvement for the Ba 
and Cu dominated modes, bringing their frequency values in excellent agreement with their 
experimental counterparts and also reproducing the mode admixture very accurately. From 
isotopic substitution experiments on the copper-sites the values for all three independent 
elements of the dynamical matrix of the Ba-Cu subsystem could be extracted.^ It can be 
observed (Table [TV| ) that the LDA calculations underestimate both diagonal terms by 15% 
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whereas the off-diagonal term lies within the experimental uncertainty. The respective 
values from GGA calculations exhibit perfect agreement with the experiment. 

With regard to the oxygen modes, we find that the frequency of the in-phase 0(2)+0(3) 
mode is increased ranging approximately 8% below the experimental value, whereas the 
phonon frequencies of the out-of-phase 0(2)-0(3) mode remains basically unaltered ac- 
curately reflecting the experimental value. The increase in the diagonal force-constant of 
the apex oxygen is not accompanied by a corresponding increase in the frequency of the 
0(4)-mode. This is due to the changes in the coupling between 0(4) and the other atoms. 
Concerning the mode admixtures, no experimental data are available to fully determine the 
off-diagonal elements of the 0(2)-0(3)-0(4) subsystem. Nevertheless, site-selective isotope 
substitution experiment S0 indicate that the apex mode is largely dominated by the 0(4) 
vibration and that the hybridized 0(2) +0(3) mode has a small but noticeable admixture of 
the apex oscillation. On contrary, the 0(2)-0(3) mode is decoupled from the 0(4) motion. 
Our calculations confirm this picture. 

VI. CONCLUSIONS 

We have determined phonon frequencies and eigenvectors of the five Raman-active A lff 
modes in YB^CuaOy employing frozen-phonon calculations on the basis of the full-potential 
LAPW method. Atomic-force values were used to establish the atomic equilibrium posi- 
tions and to map out the phonon potentials. When treating the electronic exchange and 
correlations by LDA, the calculated phonon frequencies display deviations from the mea- 
sured values by approximately 10%. The only exception is the out-of-phase 0(2)-0(3) 
mode where the agreement with experiment is excellent. Going beyond LDA by invoking 
a GGA exchange-correlation functional leads to an overall improvement of the frequencies. 
The two low-lying Ba-Cu modes as well as the 0(2)-0(3) mode are in perfect agreement 
with experimental values whereas the frequency of the 0(2) +0(3) mode is enhanced but 
still 8% too small. In comparison to LDA calculations only the 0(4) frequency remains 
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unaltered ranging 10% below its measured counterpart. Regarding the mode admixtures we 
find that our GGA results for the two Ba-Cu modes very accurately reproduce the diago- 
nal and off-diagonal force-constants which have been determined by isotopic substitution 
experiments. Additionally, our calculations show that the 0(2)-0(3) mode is completely 
decoupled from the 0(4) oscillation, whereas there is a finite contribution of the 0(4) vibra- 
tion to the 0(2) +0(3) mode. The domination of the 0(4) displacement in the apex mode 
is also revealed. These findings are compatible with Raman measurements performed on 
site-selective substituted samples. In summary, these results indicate that first-principles 
calculations based on density functional theory are an appropriate tool for the investigation 
of structural and vibronic properties of HTC cuprates, where the treatment of exchange and 
correlation effects is an important issue for this task. 
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TABLES 



TABLE I. Force contributions in mRy for a mixed distortion (-0(2), -0(4)). 





Ba 


Cu(2) 


0(2) 


0(3) 


0(4) 


position [c] 


0.1815 


0.3530 


0.3740 


0.3787 


0.1540 


f hf 


46.66 


123.82 


0.76 


-8.77 


290.52 


-pcore 


-36.08 


-88.30 


6.45 


0.75 


-188.87 


■pIBS 


-13.03 


-35.22 


6.50 


7.91 


-75.75 


pat 


-2.45 


0.30 


13.71 


-0.11 


25.90 



TABLE II. LDA phonon frequencies in cm and eigenvectors of the five Ai g modes in 
YBa2Cu307 at q = compared to experimental data00. 



Exp. 


Theory 


Rel. Deviation 


Ba 


Cu(2) 


0(2) 


0(3) 


0(4) 


116,118 


103 


-11%,-13% 


0.85 


0.52 


0.05 


0.05 


0.00 


150,145 


130 


-13%,-10% 


0.53 


-0.84 


-0.08 


-0.07 


0.06 


335 


327 


- 2% 


0.00 


0.02 


-0.81 


0.59 


-0.05 


440 


387 


-12% 


0.02 


0.09 


-0.51 


-0.74 


-0.43 


500 


452 


- 9% 


0.03 


-0.11 


0.28 


0.32 


-0.90 
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TABLE III. GGA phonon frequencies in 


cm 1 and 


eigenvectors 


of the 


five A lg 


modes in 


YBa 2 Cu 3 07 


at q = compared to experimental data00. 










Exp. 


Theory Rel. Deviation 


Ba 


Cu(2) 


0(2) 


0(3) 


0(4) 


116,118 


115 -l%,-3% 


0.91 


0.40 


0.05 


0.02 


0.02 


150,145 


144 -4%,-l% 


0.41 


-0.91 


-0.06 


-0.04 


0.03 


335 


328 - 2% 


0.01 


0.04 


-0.81 


0.60 


-0.02 


440 


405 - 8% 


0.02 


0.07 


-0.56 


-0.77 


-0.30 


500 


452 - 9% 


0.02 


-0.04 


0.20 


0.23 


-0.95 



TABLE IV. Force constants in N/m for the sublattice vibrations of Cu and Ba compared to 
measured values0. kn and k 22 denote the diagonal force constants, ki2 the coupling term. 

LDA GGA Exp. 

ku 101 119 118.6 (-5.3/+5.0) 

k 22 69 77 81.3 (-2.2/+2.5) 

k i2 -17 -15 -16.3 (-3.3/+S.6) 
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